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Abstract 

Three dimensional electromagnetic gyrofluid simulations of the ideal ballooning mode blowout 
scenario for tokamak edge localized modes (ELMs) are presented. Special emphasis is placed on 
energetic diagnosis, examining changes in the growth rate in the linear, overshoot, and decay 
phases. The saturation process is energy transfer to self generated edge turbulence which exhibits 
an ion temperature gradient (ITG) mode structure. Convergence in the decay phase is found 
only if the spectrum reaches the ion gyroradius. The equilibrium is a self consistent background 
whose evolution is taken into account. Approximately two thirds of the total energy in the edge 
layer is liberated in the blowout. Parameter dependence with respect to plasma pressure and the 
ion gyroradius is studied. Despite the violent nature of the short-lived process, the transition to 
nonlinearity is very similar to that found in generic tokamak edge turbulence. 

PACS numbers: 
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I. INTRODUCTION 



Edge localised modes (ELMs) are bursty, quasi-periodic expulsions of energy and particles 
from the plasma edge in the high-confinement state of toroidal magnetised plasmas l|-|3|. 
ELMs originate from the steep pressure gradient pedestal in the outer closed field line region 
of a tokamak. A substantial release of magnetic and potential energy, caused by a rapid rise 
of the rate of momentum transfer through the flux surfaces, leads to enhanced transport 
and loss of heat and particles into the open scrape-off layer (SOL) field line region. The 
filamented peak heat fluxes during ELMs on the bounding divertor plates present serious 
restrictions on the performance of future fusion experiments like ITER 



^ and necessitate 



the development of suitable control techniques 

ELMs generically have been observed in all divertor tokamaks since the initial discovery 
of the high- confinement "H-mode" state in the early 1980s [8.], while characteristics like 
frequency and intensity are found to vary widely depending on experimental conditions {q]. 
Stellarator experiments also are reported to show similar bursty quasi-periodic edge activity 



in the presence of edge transport barriers pi|. Indication that equilibrium drifts are part 
of the process has been found by experimental diagnosis of an asymmetry of the energy 
and particle fluxes on each of the divertor plates and evidence that this is sensitive to the 
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Hl3|. 



directions of the toroidal current and magnetic field 

Phenomenologically, edge localised mode events in tokamaks have been compared to solar 
flare eruptions l4, llSj, and, in the cataclysmic variability of observed radiation emanating 
from the plasma, they may seem to bear also some remote resemblance to the outbursts of 



pressure-driven dwarf novae 



16|, although the specific instability mechanisms behind these 



phenomena are clearly of a substantially different nature. 

On the other hand, the nature of ELMs has also not yet been completely clarified, in the 
sense that there is no first-principles based predictive theory available. A number of charac- 
teristically different variations of ELMs have been observed in experiments (e.g. J?ype I, II 
or III ELMs) which may actually be caused by different instability mechanisms [itI . "Type 
III" ELM observations show a remarkable similarity to global bursts found in computations 
of drift-wave turbulence, where the turbulence generated flows and fluxes are closely tied to 
the self-consistent equilibrium evolution in the plasma edge 
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19|. 



"Type I" ELMs in the experiment specifically show stronger magnetic activity during the 



2 



burst and are usually associated with the onset of a magneto-hydrodynamic (MHD) ideal 
ballooning mode (IBM) instability when a threshold pressure gradient is reached in the steep 
edge pedestal of an H-mode plasma [l^]. The possible role of an additional current driven 
instability ( "peeling-ballooning mode" ) for "Type I" ELMs has also received recent interest 
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22|. 



The ideal ballooning instability with nonlinear phases or aftermath has provided a 
paradigm for recent and current study of the ELM phenomenon theoretically. Explosive 
MHD instability and associated critical phenomena have been advanced analytically fl^ . 
Early Braginskii fluid simulations were given in support of this 23|]. Large, nonlinear MHD 
codes have been studying this in more detail, resolving several numerical problems and pay- 
ing specific attention to the details of the magnetic geometry, with mixed results on the 



phenomenology of the instability beyond its linear phase |24j-|28|. Whether or not the ideal 
ballooning instability or its peeling/ballooning variant is actually responsible for ELMs oc- 
curing in experiments (here noting the absence of a demonstrably well-resolved L-to-H mode 
transition in computations from first principles), the nonlinear phase of this "ideal MHD 
blowout" or "IBM blowout" phenomenon is of physical interest. The ability of the insta- 
bility to destroy the edge pedestal layer is undisputed, and the term "blowout" is germane. 
Typical physical parameters lead to a situation wherein the ideal MHD interchange growth 
rate is comparable to the parallel Alfven transit time, one and the same with the basic ideal 
ballooning criterion. This is fast with respect to MHD but not to microturbulence. Provided 
they can simultaneously treat global MHD (i.e., the spectrum covers both global and ion 
gyroradius scales), edge-turbulence computations can also address the phenomenon and are 
even in a position to treat parts of it which fall outside the paradigm of the MHD model. 

Here we present nonlinear gyrofluid computations of IBM blowout events localised to 
the edge/SOL region of a tokamak plasma, resolved to below the ion gyroradius pj. As in 
the large MHD models the self-consistent evolution of the equilibrium is not only included 
but is also an integral part of the dependent variables. An electromagnetic 6-moment gyro- 
fluid model ("GEMR") for both electrons and ions is used, including an energy- conserving 
treatment of finite ion Larmor radius effects and higher-moment terms which occur in the 
toroidal drift 29|. The energetic consistency is a key feature since the reaction of parts of 



the system such as flows and currents which have low energy content can play a central role 
in the indirect nature of nonlinear dynamics by serving as transfer channels. Although the 
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gyrofluid model has a different formulation of nonlinear polarisation (the generalisation of 
"vorticity" to a two-fluid setting) than a fluid or MHD model, full correspondence in the 
regime of validity of the latter models has been shown 30|. In this sense the model is a 
superset of global reduced MHD and electromagnetic microturbulence. The global dynamics 
is in the reduced MHD regime [3l| due to low absolute edge beta values (/3 = Svrp/i?^ ~ 10^'^) 
but the small scales depend on treating the ion gyroradius to arbitrary order. The model 
does not assume instabilities to occur at any particular scale, but when they occur the 
resulting spread of the spectrum is found to reach pi within a few eddy turnover times of 
the onset of nonlinearity in every case computed. 

The work reported herein addressed the following points. The relationship of the MHD 
dynamics to the ideal threshold is much discussed but in the gyrofluid model several mech- 
anisms of microinstability are also present. The threshold issue is obscured by the existence 
of ion temperature gradient (ITG) drift instabilities at all beta values, with a smooth tran- 
sition between them. Nevertheless, varying character in the saturation phases is observed at 
different values. We also address the issue of possible "rho-star" dependence (here the local 
rho-star is given by Ps/Lx where ps is pi evaluated at Tj = Te and Lt is the scale length 
of the Te profile, both evaluated at the mid-pedestal location in minor radius) which would 
indicate a role for diamagnetic drifts. The gradient is sufficiently steep that Ps/ Lt > Lx/qR 
so that none of these effects can be ordered small. We examine the necessary resolution to 
obtain converged cases, which gives the spectral range actively involved in nonlinear satura- 
tion. The main instability at experimentally relevant beta is the ideal ballooning instability 
but the main saturation process is energy exchange with broadband electromagnetic drift 
wave/ITG turbulence which the instability itself generates. Due mostly to this, the results 
turn out to be outside of both the magneto-hydrodynamic (MHD) and collisional Braginskii 
paradigms, on which most previous approaches have been based. On the other hand, the 
GEMR model still uses delta-f equations, so the actual profile phenomenology in the SOL 
region is not well represented, and as in other work with GEMR we concentrate on the 
properties of the nonlinear dynamics 32 1. 

Following sections give the details of the GEMR model, a discussion of the physical 
difference between ITG and MHD dynamics and their relative roles, the methods used to 
pre-set the axisymmetric equilibrium state so that the computations are done in the absence 
of axisymmetric oscillations whose decay times can be slower than the IBM rise time, the 



4 



computations of the blowout phase itself, and the possibility and limitations of quantitative 
comparison to experiment. 



II. THE GYROFLUID MODEL WITH RADIALLY DEPENDENT AXISYMMET- 
RIC GEOMETRY 



The gyrofluid model used herein is given in Ref . |[29|] . It is based on the original derivations 



given in Refs. 
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34| . with corrections given as motivated by free energy conservation. An 



alternative derivation using the conservation laws of the underlying gyrokinetic model is 
given in Ref. 35|. The equations are normalised to spatial scale a and time scale a/cg where 
a is the minor radius and Cg is the reference sound speed given by = T^jMi^ that is, in 
terms of the electron temperature (in energy units) and ion mass. The main parameters are 
the drift parameter, electron dynamical beta, and electron collisionality, respectively given 
by 

Pe = —p^ Ve = — (1) 



ps 

p* = — 
a 



where c and e are the speed of 
and Braginskii collision time 



52 CsTe 

ight and fundamental charge, Pe and are the pressure 
I of the electrons, and ps = CsMic/eB is the drift scale. 
Note that if Tj = Tg then ps = pi. The reference for Mj is the deuterium mass M^. For 
the electrons the physical value of the mass nie/Mij = 1/3670 is always used. With global 
geometry there is no single magnetic field line connection length, but one with a profile with 
2TTqRo a function of the minor radius coordinate, where Rq is a constant giving the reference 
value of the major radius. 



A. Species constants and gyrofluid moment equations 

Each species z has its own set of gyrofluid moment variable equations, one each for density 
riz, parallel velocity Uz\\, parallel and perpendicular temperature T^y and Tz±, and parallel 
and perpendicular energy components of the parallel heat flux q^w and qz±. They are coupled 
to the electrostatic and parallel magnetic potentials and A\\ through self consistent field 
equations. The tilde denotes a dependent variable to distinguish from constant parameters. 



The correspondence of these equations to the Braginskii fluid model is given in Ref. 30 |. 



Each species is characterised by a background charge density, temperature/charge ratio, and 
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mass/charge ratio, given by the normahsed parameters 

""^"^ '''^ZT, ^'^JM- 

For electrons = = — 1 and /i^ = —me/ Mi. For the main ions = /i^ = 1 while = Xj 
is kept as a parameter. For the trace ions = is always taken. For hot trace ions is 
large while ^z is still moderate. The species gyroradius pz is given by {pz/PsY = PzTz, which 
is always small for electrons, close to unity for main ions, and can be moderate or large for 
trace ions. For this study a simplified geometry is taken with B = 1 except in curvature 
terms, so magnetic pumping and shaping effects |37[ are neglected. Since dynamics at the 



electron gyroradius scale are neglected the gyroaveraging effects on A\\ are not treated. 
The moment equations are given by 

+ Vle-^Uz + W£;-VT2_L + V|jM^|| 



dt 



^(0c + ^ + r.^^4^| (3) 



^ (^/3e^|| + PzUz\^ + Pz^E-^Uz\\ + Pz^E-^qz 
+ V|| {4>G + TzPz 



^dTzW 1 ~ ^ , 

2^ + 2"i^-VT,||+V||K|| + g,||) 

= ^( 4>o + rzPz, + 2rzfz^^ (5) 



dt 



+ me-VTz^ + Wij- V(n^ + 2T2_l) + Vyg^x 

= 2 + ^ 

+ 5a (6) 



= (7) 



dt 



+ yU^Wij- V(m2|| + 2g^x) + V|| {vLg + r^.f^^ 



(8) 



The linearised pressures are given by 



Pz 



riz + T, 



Pz 



riz + T^ 



(9) 



CoUisional dissipation is controlled by a collision parameter z/^ analogous to z/g and a set of 
numerical constants for each species, 



a. 



(10) 



which are the thermal force, thermal conductivity, and viscosity coefficients involved in 



parallel dissipation, with values given by the Braginskii model 36|, though that regime is 
never reached in core turbulence. For electrons these coefficients are 0.71 and 3.2 and 0.73, 
and for ions they are and 3.9 and 0.96, respectively. The resistive dissipation is given by 



with parallel current given by 



O.SlJii + 



0.71 
~3J 



qe\\ + qe± + 0.71J| 



■^11 = y^,(^zuz\ 



(11) 



(12) 



noting that it is the electron heat fluxes that enter R^i and as with Pe^w the Rei term enters 
the parallel motion in the same way for every species. The anisotropy dissipation is given 
by 



37r^ V " 

The thermal conduction components are given by 



5/2 

5/2 ~ 
K±_ = ^izTz^z [qz± + ^Aajn 

Hz V 

5/2 _ _ 
Ka = l.^^HzTz^z {qz\\ - l-5g^±) 

Hz 



(13) 

(14) 

(15) 
(16) 
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The Landau damping effects are modeled by 

= (1 - Q.l2^q^RlV\) (17) 
qtio 

witli tliermal velocity Vz given by = Tz/pz- 

The gyroaveraging is done through Pade approximants to operators which would act in 
Fourier space, 

ri^fl-kvlV and r.^pl^ (18) 



V 2^^ V ^ "^'dpl 

and then the gyroaveraged potentials are 

0G = ri0 and l^G = r20 (19) 

noting that they are species dependent. The main ExB advection and the FLR correction 
pieces are given by 

u^;- V/ = /] - iy±^lf + i/yVj/ (20) 
Wi,-V/=[QG,/] (21) 

where the [, ] symbols denote the nonlinear bracket defined below, and iy±,iy\\ denote the 
artificial dissipation coefficients. The parallel derivative is given by 

The curvature operator is given by 

}C{f) = -[\ogB'j] (23) 



B. Field equations for polarisation and induction 

The species are coupled through two equations which give the self consistent response 
of the field potentials to the evolution of the gyrofiuid moment variables. The electrostatic 
potential is governed by quasineutrality, which sets the total space charge density to zero. 
The density for each species is given by a part due to the gyrocenters and another due to 
polarisation, which is what sets 0. This polarisation equation is given by 

= J2az tin, + r^nA (24) 
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where 

To = {i-pyir' (25) 

gives the gyroscreening of the main ions. In this version of the model the electron and trace 
ion contributions to gyroscreening are neglected, due to small mg and zero a^, respectively. 
On the right hand side the electrons and main ions enter with oppositely signed values 
and again here the trace ions do not enter. This establishes them as a trace species. 

In a similar way the trace ions are left out of the induction equation due to the zero a^, 
so that 

- p^Vil|| = J|| and J|| = ^a,u^\\ (26) 

z 

determines A\\ noting the way that the normalisation scales it with respect to p^,. 

More detail, including energy conservation and the relationship of these equations to it. 



is given in Refs. 29|, l35 |. 



C. Representation of the self consistent magnetic geometry 

In these expressions the operators and d/ds and the nonlinear brackets are determined 
by the representation of the geometry that is used. 

The magnetic geometry is a simplified representation of an axisymmetric magnetic field 
using field aligned Hamada coordinates. The field representation is 

B = /Vv5 + V^xVy? (27) 

where cp is the physical toroidal angle about the symmetry axis, is the magnetic flux 
function, and / = /(^E') is given by the constraint of MHD equilibrium. Without loss of 
generality one may define magnetic flux coordinates {V, 6, (} with the following properties, 

B^ = B^ = x,v B^ = ^y (28) 

where V = V^(^') is the volume enclosed by the surface with flux \E', then x = x(^) ^iid ip = 
4>(y) are two further flux functions found, and the subscript denotes the partial derivative 
with respect to V. The ratio q = dip/dx is another flux function which gives the pitch of 
the magnetic field. The angles are cyclic on [0, 1] which determines the coordinate Jacobian 
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to be unity. Hence the metric tensor gij satisfies det gij = g = 1. The function x is found 
by the constraint that 6 is cychc on the unit interval 

Xy'=iJ^ and e = xyfj^ (29) 



B ■ Vr/ ' J B-Vr] 

where is a simple cyclic coordinate on the flux surface in the poloidal plane (e.g., path 
length on the curve, or the physical angle position about the magnetic axis with respect to 
any branch cut). Then q is found 

^ ^ (30) 



' 27rx,v\R\ 

where the angle brackets give the flux surface average defined as 

^ \ I de I dr] 

Now, ipy = qx,v is also defined and the toroidal angle coordinate is given by 

which completes the prescription. 

The field aligning is a one-to-one and onto coordinate transform from {V, 6, to {x, y, s} 
given by 

x = V/a^ Vk = 16 — C — ak{x) s = 9 (33) 

where a is a reference minor radius and is an arbitrary function of x which is chosen 
to make the off diagonal perpendicular metric element g^^ vanish at a particular location. 
This is called the shifted metric procedure 3^. The point is that at any position in s where 
perpendicular drifts or V^^ is evaluated, the coordinate elements are rectangular, and the 
information contained in magnetic shear enters only in derivatives in the parallel coordinate. 
The magnetic field components now satisfy 

B^ = Bl = B^ = x,v (34) 

at all locations in space for any choice of a^- Hence although the perpendicular coordinate 
elements are rectangular only for one particular location in s the field aligning is exactly 
satisfied everywhere. 
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The coordinate metric elements g'^^ — Vx' • Vx^ are found in the {V, 9, (} representation. 
Then the grid locations s — Sk are all given their own members of the family of these field 
aligned coordinates via 



ak{V)^qSk+ [ dV 
Jo 



gVV 



(35) 



e=sk 



Then we have 



gxx ^ a-^g^^ (36) 

gr = (37) 
gvy^ q2gee_^qg0c + gCC (38) 

g^^= a-^g^^ (39) 

(40) 

= g'' (41) 
at s — Sk- The drift tensor elements are given by 

F = e-B with F,y = x,v (42) 

that is, Fxy is the only independent, nonvanishing component. We have 

c cF 
Vf-—BxV<P = V<P- — -Vf (43) 

which defines the bracket [(f), f] for any scalar field quantities (j) and /. With the electro- 
static potential this gives the ExB advection term, for electrons for which FLR corrections 
are neglected by taking pe — > 0. For ions this is generahsed by the FLR corrections as given 
above. 

For local equations as those used here (nonlinearities kept only in advection and magnetic 
flutter) the geometry model must be modified to retain free energy conservation by the 
equations. In particular, any function of x such as a radially varying normalised parameter 
(such as the temperature/charge ration r^) multiplying any of the curvature terms will 
cause them to fail to conserve free energy. The compressibility is already split away from 
advection and placed into the curvature terms, so the brackets must retain the properties of 
incompressible advection. This means that the model for F which is used must satisfy 

cFn 

F ^ Fo with V • = (44) 
11 



and to avoid any confusion the curvature terms are written with logS^, which can be 
given any spatial dependence. Correspondingly, the MHD version of the continuity equation 
(neglecting diamagnetic fluxes) has already been manipulated according to 

— + V ■ nu = (45) 

at 

d \ 

+ u- V logra + V- u = (46) 



and then terms such as Myb ■ V are neglected due to the ordering. The result is the same as 
the local continuity equation 

d \ 

— + u ■ V n + noV ■ u = (47) 



dt 

where no is a normalising constant, if we identify the dependent variable n/riQ with logn. It 
is important to note that with the self consistent background being evolved this also includes 
the profile: Vlogn — )■ Vn. 

The simplified geometry is now defined by its operators. The Laplacian and gyroaverg- 
ing operators neglect d/ds. The form in the Ampere's law (which does not involve the 
gyroradius) is given by 

The form in the gyroaverging operations and polarisation is given by 

dx B"^ dx B"^ dyl 



n = l'~ + %lj (49) 



The brackets are defined by 

[/, g] = Qlf, 9]xy + [/, 9]xs - q,x{s - Sk) [/, g\ys (50) 
for any two scalar fields {/,(?} with each bracket piece defined by 



^ dx^ dx^ dx^ dx^ ^ 
with the third piece not used since these are always evaluated at s = 



III. MODE CHARACTER: ITG VS. MHD 



Nonlinear computations on ideal or peeling ballooning mode ELM scenarios in tokamaks 
have previously been based on magnetohydrodynamic (MHD) models 22, |26|-|28|, and the 
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respective codes had largely been tailored towards mode structures expected by linear anal- 
ysis. Restriction to single-fluid MHD equations including two-fluid correction terms allows 
computations of large to mesoscale dynamics in realistically shaped 2D tokamak equilibria 
including X-point geometry crossing the separatrix, due to the lack of a need to use fleld- 
aligned coordinates. Recent MHD simulations are able to reproduce spatio-temporal mode 
structures, divertor footprints and other characteristics with good agreement to experimen- 
tal diagnostics, in particular for the growth and immediate blowout phase of the flnger-like 
instability 40 1. 

However, as soon as the dynamics becomes nonlinear, the spectrum broadens to include 
scales normally associated with drift wave dynamics. This turbulence is basically of a drift- 
Alfven type with strong edge ITG character (see Ref. 39j for the signatures of the various 
mode structure types), made more powerful than otherwise by the energetic access to the 
long- wavelength MHD component. The ITG character results from the steeper logarithmic 
gradients in both temperature proflles, and the lack of involvement of ion temperature 
fluctuations in the parallel nonadiabatic response of the electrons. Although edge turbulence 
of this type (cf. [19]) is strongly suppressed in the H-mode phase and is not initially involved 
in the instability phase, the experience of edge turbulence becomes relevant as these physical 
components become involved in all the nonlinear phases of the blowout. The saturation and 
aftermath of the blowout, which in an experiment carries a large part of the actual transport 
losses that lead to a degradation of the pedestal, should be expected to involve physics not 
contained in the MHD model. 

Ion temperature gradient (ITG) driven modes and ideal ballooning modes are both essen- 
tially caused by the gradient and curvature driven interchange instability and show similar 
character in their initial linear growth phases. Instability is achieved above their respective 
critical gradients, determined by the ratio 77, = Ln/Lxi > rj^ between density gradient to ion 
temperature gradient scale lengths for ITG, and by the ideal MHD ballooning parameter 
ttAf = qKVP > s for IBM, where s is the magnetic shear parameter and rjc is a critical 
threshold which depends on beta, collisionality, and toroidicity. The principal difference 
between these modes is the relative role of the parallel Alfven responses which tend towards 
enforcement of an adiabatic response in the electrons (parallel force balance, with the elec- 
tron pressure gradient). This is nonexistent in an MHD model and subdominant generally 
for an MHD instabililty. But it constrains the electrons if the aM parameter is below thresh- 
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old, still allowing the ITG instability because the adiabatic response does not involve the 
ions. The ITG instability exists at all um values, and furthermore for the edge situation the 
ITG instability and mode structure transitions smoothly to a drift wave one in a nonlinear 
setting for rji ~ 1. Hence there is actually no finite threshold in values of rji or um for the 
nonlinear edge situation. Therefore, the existence of a threshold in the experiment is the 
same as the mechanism which maintains the H-mode, which is not yet well understood. It 
follows that a gyrofiuid model computation of an ELM scenario with enough resolution to 
allow for ion-gyroradius based dynamics has to face the lack of a simple linear threshold in 
the parameter space. 

The transition from initial (micro-)instabilities to generic edge turbulence was studied in 



detail in Ref. 



4l|. As the most unstable linear modes crystallise out of an initial random 



bath of small-amplitude perturbations, the linear growth rate rises and becomes steady. The 
maximum value of the instantaneous growth rate of total fluctuation free energy E given by 
Tit) = {l/2E){dE / dt) may be taken as the maximal linear growth rate. The curve of r(t) 
then falls very sharply to zero (over about 10L_l/cs) as saturation occurs. There is some 
structural adjustment over the next few 100L_l/cs as the spectrum fllls out, and then the 
turbulence is fully developed. But over the adjustment phase the value of F is well below its 
previous maximum. We will use the same diagnostic herein, except that the total ion ExB 
heat flux Qi averaged over the computational domain is used instead of E because most of 
E is represented by the self consistent proflles while Qi is entirely due to fluctuations. 

The IBM ELM blowout scenario is similar to this, initially, except that the instability is 
not a microinstability. Nevertheless, the scale differs by less than an order of magnitude: 
the toroidal mode numbers for linear ITG instabilities are in the range of n = 30 — 100 while 
the main ideal ballooning mode is near mode number n = 10, on the entire flux surface, for 
these typical Ps/L^ values (recall ps = Pi for Tg = Tj, and the ratio Ti/Tg is not far from 
unity in experimental cases). The ITG mode numbers are determined by the ion gyroradius 
and the dimensionless scale ratio Lt / qR, while the MHD values are determined by the width 
of the pedestal. Since Lt is not larger than about 30pi in an H-mode pedestal, these scales 
are not disparate. 

Due to the closeness of the native scales, the IBM instability very quickly transfers energy 
to the ion gyroradius scale, which is only about two or three steps away in a turbulence cas- 
cade which acts at a factor of two at each step. This brings the ITG and drift Alfven physics 
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FIG. 1: Visualisation of the tokamak ELM burst: the edge region of a circular poloidal tokamak 
cross-section is shown on the left with computational perpendicular s-domains (s = 2, 5, 8, 11, 14 
of Us = 16 total) mapped onto the circle during the maximum growth phase (t = 38 in units of 
a/cs) of an ideal ballooning mode. On the right, the outer midplane area (s = 8 with size rix = 48 
and Uy = 512) is magnified and the radial domain stretched by a factor of 2 for better visibility at 
times t = 38 (max. flux) and t = 42 (turbulent aftermath). 

into play, and in any nonlinear stages the latter involves stable shear Alfven component with 
substantial energy content. The moment of saturation is defined as the time at which the ion 
ExB heat flux averaged over the computational domain is maximum; we refer to this as the 
"peak-fiux" time. The results shown below indicate that coincides with the establishment 
of the fully turbulent regime. Representation of this phase requires resolution of all scales 
and inclusion of the appropriate drift wave physics in the model. How much resolution is 
actually needed depends on the strength of the blowout, which is situation dependent. We 
therefore include a resolution scan on the toroidal mode number spectrum. The result is 
that convergence is reached in the temporal phases just after peak flux only if the spectrum 
reaches down to kyPi = 1. 

The self-consistent equilibrium coupling in our gyrofiuid model only allows treatment of 
shifted circular s — a geometry (the difference to MHD codes is that the high-resolution FLR- 
gyrofiuid dynamics necessitates field-aligned coordinates). Effects of flux surface shaping 
are postponed to later studies. Our focus here is on basic physics issues regarding the 
ability of the model to capture MHD phenomena (see also Ref. ISj concerning global Alfven 
oscillations with self consistent profiles), resolution, turbulent character, scalings and tests 
which are also accessible by s — a geometry. In particular, the theory of explosive instability 
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has been formulated for circular geometry [ij] and is therefore accessible by the present 
model. 



IV. MODELLING OF INITIAL PROFILES 

First-principles based local drift wave edge turbulence simulations are not able to obtain 
a realistic H-mode edge state with the known experimental characteristics: correct den- 
sity and temperature pedestal profiles shapes or strength of flow shear are not obtained 
by self-consistent evolution by specifying core sources only, nor has a threshold transition 
character been found in any verified edge turbulence simulation |l9|. Therefore some kind 
of "modelling" has to take place when the IBM instability (as an H-mode phenomenon) and 
its subsequent nonlinear evolution is simulated with a nonlinear gyrofluid turbulence code: 
Although the realistic development of an edge transport barrier (and thus a full ELM cycle) 
can not be directly obtained, one still may prescribe the H-mode pedestal profile before the 
onset of an ELM, known from experimental data, as an initial state for the simulation. As 
a base case for the prescribed pedestal profiles the well diagnosed edge characteristics of 
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ASDEX Upgrade H-mode shot #17151 is used here 

The local parameters, taken as mid pedestal values, correspond to electron and ion tem- 
peratures Tg = 300 eV, Tj = 360 eV, densities = Ui = 2.5 ■ 10^^ m~^, magnetic field 
strength B = 2.0 T, major torus radius R = 1.65 m, aspect ratio R/a = 3.3, perpendicular 
temperature gradient length = L± = 3.0 cm, density gradient length L„ = 6.0 cm. The 
profile of the safety factor g = 1.5 + 3.5(r/a)^ is parabolic yielding local values at the last 
closed flux surface (LCFS, r/a = 1.0) of qa = 5.0 and Sa = 1.4. The radial domain of the 
simulations covers a range of Lj_ on either side of the LCFS. The nominal pressure values 
are Pe = rifTf. and Pi = UiTi. All times are given in units of a/cg. 

The computational grid is given in terms of the {x, y, s} coordinate domains. The spacing 
is equidistant. For each grid point in the third coordinate s = Sk the y-coordinate is i/k as 
defined in Eqs. fl33|35l) so that for each grid point Sk the coordinate system y = is the 
one which has g^^ = at s = s^. For the nominal case = 512 perpendicular and Ug = 32 
parallel mesh points are used. The radial domain with = 64 spans the plasma edge 
region between the H-mode pedestal top, with plasma core parameters as inner boundary 
values, and the outer bounded scrape-off layer region (r/a = 1 ± 0.06). This represents a 
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FIG. 2: Radial tokamak edge pedestal profiles of the plasma pressure before and after an ELM 
blowout event across the separatrix into the scrape-off layer (SOL). The area shaded in grey shows 
the present radial numerical simulation domain. 



spatial range from the global scale to smaller than the ion gyroradius scale, with a ratio 
S = ps/cL = 0.001875 between ion gyroradius ps and minor torus radius a. In the SOL 
region (r > a) the parallel boundary condition is replaced by a Debye sheath model whose 



treatment is given in Ref. 



44| 



V. PROFILE PRE-EQUILIBRATION 

The initial conditions are thus prescribed and are based on experimentally diagnosed 
radial temperature and density pedestal profiles T(r) and n(r) for each species (electrons 
and ions). A consistent electrostatic potential 0(r) is derived by numerically solving the 
neoclassical balances (parallel dynamics, toroidal drifts, and coUisional dissipation, but not 
transport) in a pre-processing step with a modified (zonally frozen) GEMR setup, resulting in 
a time-steady 2D dissipative solution. The s — a MHD equilibrium in turn is set internally in 
GEMR by computing the Shafranov shift, and in the present setup would not be consistently 



described by prescribed in terms of external or coupled shaped equilibrium solvers [18 1. 

The parallel and perpendicular electron and ion temperatures, Te\\, Te±, Tq and Ti±, 
are directly adopted and fixed from experimentally derived values by filtering the zonal 
component out of the total time derivative dT/dt = St- Their values are zonally frozen by 
setting 

dT 

-^ = St- (St) (52) 
where {St) is the zonal (parallel and perpendicular) average of all right-hand-side terms St 
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in the temperature evolution (Eqs. 5,6) including numerical dissipation. 

On the other hand, gyrocenter densities Ue and rii have to be set to obey relaxation 
relations that allow the vorticity to freely evolve into equilibrium. This is achieved by 
freezing the zonal component of the sum 11^ + r^nj, where = Ti/T^, as the density part 
of the pressure during the equilibration phase, but allowing the difference (i.e., vorticity) to 
evolve freely. This ensures that the contribution of densities to the total pressure is zonally 
frozen through the relation 

d 

— (rie + nrii) ^ Se + nSi - {Se + nSi) (53) 
while the densities relax regarding to 

^ = Se--^{Se + nS,) (54) 

at 1 + Tj 

and 

^ = S,--^{Se + nS,) (55) 

Ot l + Ti 

The numerical solution of the equilibration phase, starting directly from realistically steep 
pedestal profiles To(r), into steady state is delayed by long, weakly damped global geodesic 
Alfven oscillations. Convergence is expedited by ramping up all of the gradients gradually 
from zero to prescribed value over the first At = Tr = 50 a/cg of the run by 

^^ST-{ST) + yTo{x). (56) 

This pre-processing equilibration phase is run until convergence with reduced perpen- 
dicular resolution {nx,ny,nz)—{64: x 4 x 32), and without the ExB and magnetic flutter 
nonlinearities, which allows establishment of the 2D structure in a smooth manner. Then, 
the resolution is increased to the nominal values (64 x 512 x 32), and a random turbulent 
bath with relative amplitude 10~^ Ps/L± is added to the background pedestal profiles in- 
side the closed flux surface region. This procedure reduces transient Alfvenic and geodesic 
acoustic ringing and prepares a reproducible initial state. However, the following sudden 
release of the nonlinearities also leads to transient oscillations. Depending on parameters, 
these may still be present at the onset of the instabilities, and can obscure a clean view on 
the nonlinear growth rates, which will be relevant for the discussion below on diagnosing 
linear or explosive instability. 
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A schematic sketch of the equihbrium pedestal profiles, representing an idealised ASDEX 
Upgrade H-mode scenario 42] , is shown in Fig. [2] together with the final relaxed state after 
the ELM blowout phase (which is discussed in the next section). 



VI. COMPUTATION OF THE ELM BLOWOUT 

When this initialised pedestal pressure profile is ideal ballooning unstable, the IBM in- 
stability in GEMR simulations is observed to be linearly growing in the pedestal region and 
at the onset of nonlinearity further overshoots and saturates, representing the beginning of 
the turbulent blowout phase during which a substantial fraction of the pedestal energy is 
thrown onto the SOL. Subsequent evolution involves turbulence in both the pedestal and 
SOL regions as the original pedestal energy is dissipated. 

Previous nonlinear approaches on ELM ideal ballooning mode burst computations have 
treated only the initial growth and nonlinear phases, focused on low-wavelength modes and 
resolved only the MHD-relevant scales, excluding treatment or discussion of the ion gyro- 
radius scales. This however precludes the development of fully developed microturbulence 
which is caused by the onset of nonlinearity (robust transfer of free energy to smaller scales) . 
Ultimately, this microturbulence decides both the transition to nonlinearity and the eventual 
saturation; that is, most of the post-peak transport curve. This affects the MHD scenario 
of nonlinear explosive growth which most previous work advances. However, we find that 
the exclusion of the ion gyroradius scales produces an under-resolved situation. 

We measure the quantitative character of the growth in both linear and nonlinear phases 
with a growth curve 

defined in terms of the heat flux Qi instead of the fluctuation free energy as explained above. 
The heat flux is computed as a zonal (flux surface) average, 

Qi{x) = ^ dyds {O.bpq + Pii_)uE + {Pi± + 2nifi±)wE (58) 

^^^^ 

and then the time trace Qi{t), normalised in terms of PeCg at nominal parameters, is com- 
puted as a volume average over the central half of the radial domain. 
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FIG. 3: Total energy E{t) (top figure) and ion heat flux Qi{t) (bottom figure) during tlie ELM 
blowout: witliout source maintenance (dashed curve); with the source term switched off at the time 
of maximum instability (here: at t = 40); and with the source maintained for the whole simulation 
time. 

The IBM instability, which follows the random seeding of the pre-processed pedestal 
profile in our computations, is very violent, growing for the nominal case at a rate I — 
0.18cs/L±, just below the ideal interchange rate. On the other hand the growth curve 
r{t) appears qualitatively the same as in basic turbulence cases. At all time points in the 
nonlinear phase r{t) is well below its previous maximum. Just after initial saturation there 
is some nonlinear evolution in which r{t) crosses zero before settling down into long-term 
statistical saturation. At late times the initial blowout no longer imprints the results: with 
a fixed source one merely finds bursty turbulence thereafter. Hence, there is no evidence for 
explosive instability in this nominal case, which, if present, would be visible in a finite time 
singularity in the fluctuating free energy E ~ (5/)^ ~ (t — tcrit)~^ rising with a power p. 

The burst and the resulting decaying turbulence act to transport the pedestal plasma 
across the LCFS, where energy and particles are lost within short times by parallel boundary 
outflow to the scrape-off layer (SOL) limiter. Without maintenance of a heat and particle 
source at the core boundary, this leads to decay of the initial gradient with a rapid onset at 
the initial transition to nonlinearity at about t = 40. This time point also corresponds to 
the peak-flux time. 

However, the ion temperature profile already starts decaying during the low-amplitude 
linear growth phase, due to neoclassical (finite Larmor radius) transport (for detail on how 
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FIG. 4: Consistency check for different perpendicular resolutions by varying the grid point number 
Uy between 64 and 512. Convergence is seen for ny = 256 and higher. 



this works in a fluid model see Ref. 



43|). This situation is shown in Fig. [3] (left frame, 



dashed line), where the total pedestal energy is reduced by a quarter of its initial content 
before the peak flux phase. Therefore we test sensitivity against sources. Maintenance of the 
ion temperature profile by a fixed ion heat source localised to the inner boundary (following 
the method for driven cases in Ref. 45|]) ensures nearly unchanged gradients during the 
linear growth phase. At the end of the simulation time (t = 90) after the ELM burst, the 
energy content is reduced to around 10% of the initial content when the source is absent 
or is switched off at the peak flux time (t = 40), and to around 20% when the source is 
maintained for the whole simulation time. 

The influence of either switching off the source or maintaining the zonal profiles on the 
nonlinear ion heat flux Qi{t) is shown in Fig. [3] (right): the linear growth phase and the 
peak flux is not significantly changed by the source. As expected, the heat flux saturates 
for late times with maintained core boundary inflow and is otherwise decaying. 

Consistency checks on the perpendicular resolution are made by varying Uy between 64, 
128, 256 and 512. The evolution of Qi{t) for these cases is shown in Fig. HI The Uy = 64 
case is clearly under- resolved and overestimates the total energy decay by a factor 2. For 



128 the initial part of the post-peak decay phase disagrees with the Uy = 512 case. 



Convergence is found for ny = 256 and higher. The nonlinear saturation phase thus requires 
resolving the ion gyroradius scale. This is the clearest indication that turbulence for toroidal 
mode numbers beyond about 30 extending down to the kyPi ~ 1 range is involved in the 
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FIG. 5: Perpendicular mode number kyPs spectra for t = 38 (top) and t = 42 (bottom) around 
the peak flux phase. The vorticity spectrum (bold lines) is already flattening down to the ion 
gyroradius scale for t = 42 only a few ps after the maximum linear growth phase. 

blowout saturation process. 

The transition from linear instability to turbulence is studied during the peak flux phase 
around t = 40. Fig. |5] shows perpendicular mode number kyp^ spectra of the squared am- 
plitude of various fluctuating plasma quantities (density n, ion temperature T,, electrostatic 
potential and vorticity u) for t = 38 and t = 42. This time difference corresponds to 
At = 15/is in physical units and around 70L^/cs in local drift units, which is only slig htly 
faster than the overshoot and saturation times known from edge microinstability cases 41 |. 
Initially, the ion temperature gradient (ITG) driven microinstability and the ideal balloon- 
ing mode (IBM) compete in growth out of the random low-amplitude bath. The ITG mode 
grows strongest near the separatrix due to radially local steepening by parallel SOL diffusion, 
which in our simulations may be overestimated by using the standard fluid Bohm outflow 
boundary conditions. For t = 38 the linear IBM is clearly dominant near a toroidal mode 
number of 9-10 for the nominal AUG parameters, consistent with experimental observations 
461]. Around t = 42 rapid formation of a turbulent cascade range in the spectra is observed 
and the vorticity spectrum is already spread out to the ion gyroradius scale [kyPs = 1). 
This is a manifestation of the role of self generated drift wave turbulence in the saturation 
process. Because of the way bioth the diffusive mixing and vorticity scattering nonlinearities 
enter the drift wave physics [39i] the spectrum is held together as a unit, all scales down to pi 
are involved in the saturation phase of the overall ELM blowout transport. The involvement 
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FIG. 6: Competing growth of ITG, microtearing and IBM instabilities determine onset times and 
intensities of the blowout time traces Qi{t) when scaUng with the local plasma beta: for /3e = 10~^ 
only an ITG micro-instability is growing around t = 50 and saturating on a low H-mode transport 
level. For /3e = 2 • 10~^ the IBM instability dominates and leads to a clean ELM signal around 
t = 20. In the nominal case (/3e = 2 • 10^^) ITG and IBM MHD growth compete initially with 
similar growth rates but the ideal ballooning mode takes over around t = 30 and further determines 
the blowout. For /3e = 8 • 10~^ a microtearing instability is saturating on a gyro-Bohm transport 
level around t = 20 and is transformed into more violent ITG MHD turbulence after deterioration 
of the initial electron temperature gradient (t > 50) . 

of Pi at unit order, in turn, necessitates a nonlinear gyrofiuid model (and at some later date, 
a nonlinear gyrokinetic model). 

The characteristics of the blowout at various plasma beta is studied for 0.5 ■ 10~^ < f3e < 
8- 10~^, with the absolute /3 a factor of 4.8 higher due to the Ti/Te = 1.2 value. A competition 
between unstable growth of ITG high-n modes and IBM low-n modes is observed: For low 
P an early onset of ITG turbulence degenerates the profile and influences the strength of 
the IBM. For large /3 the nonlinear saturation is acting faster and stronger. The two effects 
lead to a rather nondistinctive result concerning the peak ELM flux as shown in Fig. [6l with 
a continuous transition between ITG and IBM triggered onset of the turbulent aftermath. 

The presence or absence of background current gradient terms, which in this model are 
set by replacing (J|| — )■ J|| + Jg) everywhere the electron parallel velocity z/|| appears in 
the equations, with Jq given by the q profile, was found to have no discernible effect on 
the result. Indeed self-consistent inclusion of VJq effects with possible role in the L-to-H 
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transition was among the motivations of the GEMR model in the first place 18|]. However, 
the magnitude of even an impressive pedestal current peak is only one to three times the 
nominal saturation current given by ngCCs, while values closer to neeCgqR/ L± are required 
to enter the energetics effectively. Since qR/L± ~ 200 for the nominal case, current gradient 
effects are very weak. 



VII. CONCLUSIONS 



The main conclusion of this study is that the qualitative nature of the saturation and 
aftermath of the initial IBM blowout is the same as for generic edge turbulence given a small- 
amplitude start. The transition from linear mode structure and energetics to turbulence 



found for these blowout cases is the same as in Ref. [4l|]. Only the nature of the linear mode 
itself differs. The blowout saturates upon its own self generated drift- Alfven turbulence, with 
a strong ion temperature component given the gradients. The vorticity spectrum reaches 
quickly to the ion gyroradius (pi) scale, requiring the gyrofiuid model and explaining the 
numerical difficulties seen with Braginskii models — it can be argued that the cases given in 
Ref. crash on entry to the nonlinear stage. Convergence in the aftermath requires 
resolving at least pi. Unfortunately, due to the lack of a self consistent H-mode state 
in a well resolved computation, no threshold is found. At lower beta values one simply 
finds generic edge turbulence driven by the temperature gradients. It is not clear that this 
scenario really describes an actual ELM, although the energetic growth and decay curves 
are not unrealistic. However, with a well resolved transition to nonlinearity in both the 
energetic peak and aftermath phases, our studies find no evidence of an explosive MHD 
phenomenology. Indeed it can be argued that nonlinear MHD processes are pre-empted 
by the efficient transfer to smaller scales through the two-fluid drift wave physics more 
commonly associated with microturbulence. 
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